PLot trial ratings
Tory Madden & Peter Kamerman
22 Feb 2018
Import data
# Import data
df <- readr::read_rds('./data/Juliane_raw_FEST.rds')Quick look
# Quick look
dim(df)## [1] 2735 11head(df)## # A tibble: 6 x 11
## # Groups: id [1]
## id block_number trial_no trial_type perc_trial_type CS US
## <dbl> <chr> <int> <chr> <chr> <chr> <chr>
## 1 8.00 baseline 1 CS_only CS_only CSplus <NA>
## 2 8.00 baseline 2 CSminus_only CSminus_only CSminus <NA>
## 3 8.00 baseline 3 CS_only CS_only CSplus <NA>
## 4 8.00 baseline 4 CS_only CS_only CSplus <NA>
## 5 8.00 baseline 5 CSminus_only CSminus_only CSminus <NA>
## 6 8.00 baseline 6 CSminus_only CSminus_only CSminus <NA>
## # ... with 4 more variables: rating <int>, rating_time <int>,
## # button_correct <chr>, button_time <int>tail(df)## # A tibble: 6 x 11
## # Groups: id [1]
## id block_number trial_no trial_type perc_trial_type CS US
## <dbl> <chr> <int> <chr> <chr> <chr> <chr>
## 1 30.0 test3 137 CSminus_USlow CSminus_USlow CSminus low
## 2 30.0 test3 138 CSminus_only CSminus_only CSminus <NA>
## 3 30.0 test3 139 CSplus_UShigh CSplus_UShigh CSplus high
## 4 30.0 test3 140 CSminus_only CSminus_only CSminus <NA>
## 5 30.0 test3 141 CSminus_USlow CSminus_USlow CSminus low
## 6 30.0 test3 142 CSplus_UShigh CSplus_UShigh CSplus high
## # ... with 4 more variables: rating <int>, rating_time <int>,
## # button_correct <chr>, button_time <int>glimpse(df)## Observations: 2,735
## Variables: 11
## $ id <dbl> 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8...
## $ block_number <chr> "baseline", "baseline", "baseline", "baseline"...
## $ trial_no <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,...
## $ trial_type <chr> "CS_only", "CSminus_only", "CS_only", "CS_only...
## $ perc_trial_type <chr> "CS_only", "CSminus_only", "CS_only", "CS_only...
## $ CS <chr> "CSplus", "CSminus", "CSplus", "CSplus", "CSmi...
## $ US <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "low",...
## $ rating <int> 0, -7, -7, -24, -32, -39, -17, -33, -40, -33, ...
## $ rating_time <int> 4824, 4342, 5833, 4173, 2689, 3275, 3460, 3201...
## $ button_correct <chr> "correct", "correct", "correct", "correct", "c...
## $ button_time <int> 8407, 5485, 7576, 5572, 3792, 4610, 4875, 4152...Plot ratings by trial type
Ratings on the FEST are plotted for each participant and then for the whole sample. First we plot by trial type delivered, then by trial type perceived.
# Make plots
ggplot(data = df) +
aes(x = trial_no,
y = rating,
colour = trial_type) +
geom_point() +
scale_colour_brewer(palette = "Set1") +
geom_hline(yintercept = 0, linetype = 2) +
geom_vline(xintercept = 10.5, linetype = 3) +
geom_vline(xintercept = 34.5, linetype = 3) +
facet_wrap(~ id) +
labs(title = "Ratings over course of procedure",
subtitle = "Faceted by ID")
df_test <- df %>% filter(trial_type == "CSplus_UStest" | trial_type == "CSminus_UStest")
ggplot(data = df_test) +
aes(x = trial_no,
y = rating,
colour = trial_type) +
geom_point() +
scale_colour_brewer(palette = "Set1") +
geom_hline(yintercept = 0, linetype = 2) +
geom_vline(xintercept = 70.5, linetype = 3) +
geom_vline(xintercept = 34.5, linetype = 3) +
geom_vline(xintercept = 105.5, linetype = 3) +
facet_wrap(~ id) +
labs(title = "Ratings over course of test phases, UStest trials only",
subtitle = "Faceted by ID")
ggplot(data = df_test) +
aes(x = as.factor(trial_type),
y = rating,
colour = trial_type) +
geom_jitter(height = 0) +
geom_boxplot() +
scale_colour_brewer(palette = "Set1") +
facet_wrap(~ id) +
labs(title = "Ratings of test stimuli according to actual trial type",
subtitle = "Faceted by ID")
ggplot(data = df_test) +
aes(x = as.factor(perc_trial_type),
y = rating,
colour = perc_trial_type) +
geom_jitter(height = 0) +
geom_boxplot() +
scale_colour_brewer(palette = "Set1") +
facet_wrap(~ id) +
labs(title = "Ratings of test stimuli according to perceived trial type",
subtitle = "Faceted by ID")
df_test %>% filter(!is.na(trial_type)) %>%
ggplot(data = .) +
aes(x = as.factor(trial_type),
y = rating,
colour = trial_type) +
geom_jitter(height = 0) +
geom_boxplot() +
scale_colour_brewer(palette = "Set1") +
labs(title = "Ratings of test stimuli",
subtitle = "According to actual trial type")
df_test %>% filter(!is.na(perc_trial_type)) %>%
ggplot(data = .) +
aes(x = as.factor(perc_trial_type),
y = rating,
colour = perc_trial_type) +
geom_jitter(height = 0) +
geom_boxplot() +
scale_colour_brewer(palette = "Set1") +
labs(title = "Ratings of test stimuli",
subtitle = "According to perceived trial type")
Prediction interval plots
Here, we plot the prediction interval for the ‘control’ stimulua type (e.g. CSminus/UStest) and depict the trials of the comparator trial type (e.g. CSplus/UStest) that fall within and outside the predictio interval. We do this for the test trials, as delivered and as perceived, and then for the acquisition trials.
# Prediction intervals for test trials as DELIVERED
df_prediction <- df_test %>%
# Filter for CSminus_UStest
filter(trial_type == 'CSminus_UStest' & id != 23) %>%
# Calculate 95% prediction interval
group_by(id) %>%
dplyr::summarise(lower_pi = mean(rating, na.rm = TRUE) - 1.96 * (sd(rating, na.rm = TRUE)),
upper_pi = mean(rating, na.rm = TRUE) + 1.96 * (sd(rating, na.rm = TRUE))) %>%
# Truncate PIs at +50 and -50
mutate(lower_pi = ifelse(lower_pi < -50,
yes = -50,
no = lower_pi),
upper_pi = ifelse(upper_pi > 50,
yes = 50,
no = upper_pi)) %>%
# Join back with df_test
left_join(df_test) %>%
# Is CSplus_UStest within the prediction interval
mutate(CC_effect = ifelse(rating > upper_pi,
yes = 'yes',
no = 'no')) %>%
# Filter for CSplus_UStest
filter(trial_type == 'CSplus_UStest')
ggplot(data = df_prediction) +
aes(x = as.factor(id),
y = rating) +
geom_jitter(aes(colour = CC_effect),
height = 0) +
geom_errorbar(aes(ymin = lower_pi,
ymax = upper_pi)) +
ylim(-50, 50) +
geom_hline(yintercept = 0, linetype = 2) +
scale_colour_brewer(palette = "Set1") +
coord_flip() +
labs(title = 'Prediction interval for CSminus_UStest trial ratings',
subtitle = 'Points are individual CSplus_UStest trial ratings',
y = 'FEST rating',
x = 'Participant ID')
# Prediction intervals for test trials as PERCEIVED
id_19 <- df_test %>% filter(id == 19)
plot(id_19$rating, id_19$trial_no)
df_prediction_perc <- df_test %>%
filter(!is.na(perc_trial_type)) %>%
# Filter for CSminus_UStest
filter(perc_trial_type == 'CSminus_UStest' & id != 23 & id != 19) %>%
# Calculate 95% prediction interval
group_by(id) %>%
summarise(lower_pi = mean(rating, na.rm = TRUE) - 1.96 * (sd(rating, na.rm = TRUE)),
upper_pi = mean(rating, na.rm = TRUE) + 1.96 * (sd(rating, na.rm = TRUE))) %>%
# Truncate PIs at +50 and -50
mutate(lower_pi = ifelse(lower_pi < -50,
yes = -50,
no = lower_pi),
upper_pi = ifelse(upper_pi > 50,
yes = 50,
no = upper_pi)) %>%
# Join back with df_test
left_join(df_test) %>%
# Is CSplus_UStest within the prediction interval
mutate(CC_effect = ifelse(rating > upper_pi,
yes = 'yes',
no = 'no')) %>%
# Filter for CSplus_UStest
filter(perc_trial_type == 'CSplus_UStest')
# Save file df_prediction_perc
write_rds(x = df_prediction_perc,
path = './data/Juliane_raw_as_perc.rds')
ggplot(data = df_prediction_perc) +
aes(x = as.factor(id),
y = rating) +
geom_jitter(aes(colour = CC_effect),
height = 0) +
geom_errorbar(aes(ymin = lower_pi,
ymax = upper_pi)) +
ylim(-50, 50) +
geom_hline(yintercept = 0, linetype = 2) +
scale_colour_brewer(palette = "Set1") +
coord_flip() +
labs(title = 'Prediction interval for ratings of trials perceived as CSminus_UStest ',
subtitle = 'Points are individual CSplus_UStest (percieved) trial ratings',
y = 'FEST rating',
x = 'Participant ID')
# Prediction intervals for all reinforced (CS+/USH or CS-/USL) trials
df_prediction_reinf <- df %>%
# Filter for CSminus_USlow
filter(trial_type == 'CSminus_USlow' & id != 23) %>%
# Calculate 95% prediction interval
group_by(id) %>%
dplyr::summarise(lower_pi = mean(rating, na.rm = TRUE) - 1.96 * (sd(rating, na.rm = TRUE)),
upper_pi = mean(rating, na.rm = TRUE) + 1.96 * (sd(rating, na.rm = TRUE))) %>%
# Truncate PIs at +50 and -50
mutate(lower_pi = ifelse(lower_pi < -50,
yes = -50,
no = lower_pi),
upper_pi = ifelse(upper_pi > 50,
yes = 50,
no = upper_pi)) %>%
# Join back with df_test
left_join(df) %>%
# Is CSplus_UShigh within the prediction interval
mutate(diff_intensities = ifelse(rating > upper_pi,
yes = 'yes',
no = 'no')) %>%
# Filter for CSplus_UShigh
filter(trial_type == 'CSplus_UShigh')
ggplot(data = df_prediction_reinf) +
aes(x = as.factor(id),
y = rating) +
geom_jitter(aes(colour = diff_intensities),
height = 0) +
geom_errorbar(aes(ymin = lower_pi,
ymax = upper_pi)) +
ylim(-50, 50) +
geom_hline(yintercept = 0, linetype = 2) +
scale_colour_brewer(palette = "Set1") +
coord_flip() +
labs(title = 'Prediction interval for CSminus_USlow trial ratings',
subtitle = 'Points are individual CSplus_UShigh trial ratings',
y = 'FEST rating',
x = 'Participant ID')
Conclusion about calibration
Calibration was poor. However, participants 11, 15 and 26 have a prediction interval for their CS-/USlow trials that does not cross zero, and they also show no CC effect according to the previous plot.
Session information
sessionInfo()## R version 3.4.3 (2017-11-30)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: OS X El Capitan 10.11.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] bindrcpp_0.2 RColorBrewer_1.1-2 magrittr_1.5
## [4] dplyr_0.7.4 purrr_0.2.4 tidyr_0.7.2
## [7] tibble_1.4.2 ggplot2_2.2.1.9000 tidyverse_1.1.1
## [10] readr_1.1.1
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.15 cellranger_1.1.0 pillar_1.1.0
## [4] compiler_3.4.3 plyr_1.8.4 bindr_0.1
## [7] forcats_0.2.0 tools_3.4.3 digest_0.6.15
## [10] lubridate_1.6.0 jsonlite_1.5 evaluate_0.10.1
## [13] nlme_3.1-131 gtable_0.2.0 lattice_0.20-35
## [16] pkgconfig_2.0.1 rlang_0.1.6.9003 psych_1.7.8
## [19] cli_1.0.0 yaml_2.1.14 parallel_3.4.3
## [22] haven_1.1.0 withr_2.1.1.9000 xml2_1.1.1
## [25] httr_1.3.1 stringr_1.2.0 knitr_1.17
## [28] hms_0.3 rprojroot_1.2 grid_3.4.3
## [31] glue_1.1.1 R6_2.2.2 readxl_1.0.0
## [34] foreign_0.8-69 rmarkdown_1.6 modelr_0.1.1
## [37] reshape2_1.4.3 backports_1.1.1 scales_0.5.0.9000
## [40] htmltools_0.3.6 rvest_0.3.2 assertthat_0.2.0
## [43] mnormt_1.5-5 colorspace_1.3-2 labeling_0.3
## [46] utf8_1.1.3 stringi_1.1.5 lazyeval_0.2.1
## [49] munsell_0.4.3 broom_0.4.2 crayon_1.3.4